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Abstract 

Solitons of a discrete nonlinear Schrodinger equation which includes the next-nearest-neighbor interactions are stud- 
ied by means of a variational approximation and numerical computations. A large family of multi-humped solutions, 
including those with a nontrivial phase structure which are a feature particular to the next-nearest-neighbor interac- 
tion model, are accurately predicted by the variational approximation. Bifurcations linking solutions with the trivial 
and nontrivial phase structures are also captured remarkably well, including a prediction of critical parameter values. 
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1. Introduction 

It has long been known that the discrete nonlinear Schrodinger (DNLS) equation is a relevant model for 
a wide range of applications including nonlinear optics (waveguide arrays), matter waves (Bose-Einstein 
condensates trapped in optical lattices) and molecular biology (modeling the DNA double strand). One of 
the reasons this model and its variants are relevant in many areas is the extensive range of phenomenology 
that the equations encompass, including the discrete diffraction, gap solitons, Pcicrls-Nabarro potentials, 
lattice chaos, Anderson localization, snaking and modulation instabilities, among other effects, see the recent 
book [1] for a review. 

In this work we consider the DNLS equation which allows for linear coupling to additional (than just 
nearest) neighbors. Such "nonlocal" interactions have been studied before, and, in particular, solutions with 
a nontrivial phase distribution [2] were identified, while a sufficiently slow decay of the interaction strength 
was found to lead to bistability of the fundamental soliton solutions [3,4]. As concerns physical realizations 
of such interactions, they are relevant in models of waveguide arrays that are aligned in a zigzag pattern 
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Table 1 

Examples of trivial (a)-(e) and nontrivial (f) phase distributions referred to throughout the text. Dashes represent non-excited 
sites, hence they do not carry phase. The corresponding solutions, which have the form of <f>„ = exp(i9„) in the anti-continuum 
limit, persist for e ^ 0. Parameters are k2 = 0.6, fci = 1.0 and fco = —2(^2 + fci). The n = ±1 phase values listed for (f) are 
approximate ones. Their exact values are given by Eq. (6). 

[5] , and in modeling the charge transport in biological molecules [6] . For other applications of the nonlocal 
DNLS systems see Ref. [2] and references therein. 



The nonlocal DNLS equation has the general form 
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where w„(i) is the complex discrete wave field, n is the integer lattice coordinate, and the real parameter e is 
the coupling strength, the coupling matrix being composed of real symmetric elements fc„m. The Hamiltonian 
and power. 
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where \\u{t, = X]?iez Wnit)\'^, are conserved quantities of Eq. (1) if the coupling matrix is symmetric. 
Solutions of Eq. (1) are obviously invariant against the phase shift, it„(t) o ?i„(t) exp(i/3) with real /3 G R 
(the gauge invariance), and against the reflection transformation, u„(t) <->■ u_„(i). The case of knm = 
Sm.n±i — 25m. n, whcrc (5m, n IS the Krouecker's delta-symbol, corresponds to nearest-neighbor interactions 
(discrete Laplacian) of the standard DNLS equation. The goal of this work is to develop a variational 
approximation (VA) for localized solutions (discrete solitons) of Eq. (1) with the extended linear coupling. 
Predictions provided by the VA will be verified by comparison with numerical results. The manuscript is 
organized as follows. In Sec. 2 wc present the equations to be solved for the computations of the steady 
states, their phase condition for existence from the anti-continuous limit, their dynamical reduction to a four- 
dimensional map, and the stability and bifurcations for the basic type of solutions that we will consider. 
In Sec. 3 we describe in detail the VA employed to describe the main type of discrete solitons supported 
by the next-nearest-neighbor coupling that include trivial phase (all nodes with same or opposite phases) 
and nontrivial phase configurations. Finally, in Sec. 4 we present our conclusions and potential avenues for 
future work. 

2. Steady states 



A natural starting point is to consider steady-state solutions, in the form of 

Unit) = 0„e'* (3) 

where amplitudes 0^ may be complex, and rescaling was used to fix the frequency as —1. Upon the substi- 
tution of expression (3) into Eq. (1), we arrive at the stationary problem 
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Throughout the manuscript, we will consider the next-nearest-ncighbor (NNN) coupling as an example to 
illustrate salient features of the nonlocal model. In addition, we assume that the coupling matrix is symmetric 
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and the lattice is uniform, hence the NNN variant of Eq. (4) is, with obviously redefined elements of the 
matrix: 

(5) 

-e(fe0ri-2 + fci0„_i + fco(/)„ + fci0„+i + k2(j)n+2)- 



2.1. The phase condition 



In the anti-continuum (AC) limit, which is defined by e = 0, solutions are defined by respective sets of 
excited sites (with a nonzero field at them), taken as 0„ = e'^" where On are arbitrary phases. It was shown 
in Rcf. [2] that the solutions initiated at the AC limit persist, as the inter-site couplings are turned on 
(e 7^ 0), if the following conditions on the phases are satisfied: 

knm sin(6'„ ^Oyyi)^^- (6) 

n^ra 

In the NNN case, Eq. (6) reduces to either 

(i) sin(6'„ — Om) = 0, for all n,m Q {1, 2, 3} or 

(ii) 6'_i —6q = 6*0 — 6*1 and cos(6'_i — 6*0) = — fci/(2fc2). In general, we consider the relative phases as trivial if 
they are integer multiples of tt, and nontrivial otherwise. In this sense, case (i) is trivial and (ii) is nontrivial. 
If a solution is composed of trivial relative phases, we say the solution is trivial (not to be mistaken for the 
zero solution, which we ignore), and likewise for the nontrivial case. 



2.2. Dynamical reduction 



Equation (4) may be rewritten as a system of first-order difference equations. For the NNN coupling, the 
fom'th order recurrence relation (5), reduces to the following four coupled first order difference equations: 

n + 1— -'n + -;—<Pn + 1" 
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Discrete soliton solutions of Eq. (5) correspond to homoclinic orbits of the fixed point at the origin, in terms 
of map (7). Decay rates of the solutions are given by A'"' , where A is the eigenvalue of the Jacobian of system 
(7) evaluated at the origin and corresponding to the stable manifold. A detailed description of dynamical 
reductions can be found in Ref. [1, Chapter 11]. For our purposes, the only information needed from this 
reduction are the decay rates. The actual computation of the manifolds is a delicate issue, since the two- 
dimensional manifolds are embedded into the four-dimensional space for real 0„, and the eight-dimensional 
space for complex <f>n ■ Standard methods for the numerical computation of the manifolds will fail in general 
due to the existence of additional, more dominant, eigen-directions. Although the detailed computation 
of the manifolds falls outside the scope of the present work, it would be interesting to approximate the 
manifolds with an appropriate parametrized cubic polynomial following a method similar to that developed 
in Ref. [7] for the nearest-neighbor (local) DNLS. 

2.3. Linear stability 

The linear stability can be analyzed in the usual way, assuming the perturbed solution as 

Unit) = + {vn + iwn)e^' + « + «Oe^'*) (8) 
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which leads to the respective spectral problem, 



We are looking for nonzero eigenvectors, i.e., solutions of the linearized system in the P{Z,C'^) space. 
The corresponding steady-state solution is called unstable if there exists at least one eigenvector for which 
Rc(A) > 0. The only solutions that are stable (for sufficiently weak coupling) in the nearest-neighbor DNLS 
equation are those with consecutive tt phase differences between adjacent sites [8,9]. As we explain below, 
the extended coupHng affects the stability (see Fig. 1). 
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Fig. 1. (Color online) Top: The power versus the coupling strength for solution of types (a)-(e) from Table 1 for the near- 
est-neighbor (local) DNLS equation, i.e. with k2 = 0. Dashed red and solid blue lines correspond to unstable and stable 
solutions, respectively. Bottom: The same branches in the NNN DNLS equation, with k2 = 0.6. Solution (f) with a nontrivial 
phase distribution (black thicker part of the line) appears in the latter case, with power which is very close to that of solution 
(e) with the trivial phase distribution. This bifurcation can be better seen by comparing the phases, see Fig. 3. The unlabeled 
(top) branch (which corresponds to initial phases of {0, 0, vr, 0, 0}) merges with the branch labeled (e) in both panels (although 
the collision occurs outside the plotting region for the bottom panel). 
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2.4. Bifurcations 

From Eq. (6) we see that in the standard nearest-neighbor DNLS equation, only solutions with trivial 
phase distributions persist for non-zero coupling. These solutions can be continued in e and gradually vanish 
through saddle-node bifurcations [10] ^ , with only the single-site (^o = 0) and two-site (^o = 6*1 = 0) solutions 
persisting toward the continuum limit (with both approaching the soliton solution of the continuous NLS 
equation). This bifurcation scenario is shown in the top panel of Fig. 1 for configurations (a)-(e) of Table 1. 
The extended coupling allows for additional types of solutions with nontrivial phase patterns, which, in 
the case of the NNN system, means one additional waveform, indicated by the thick black solid line in 
the bottom panel of Fig. 1. In these diagrams, the power is plotted versus coupling strength e. For this 
choice of parameters, the nontrivial branch (f) and the trivial one (e) lie very close to each other, making 
the identification of the bifurcation points difficult. The role of the nontrivial- phase solutions can be better 
seen through the comparison of the phases, see below Figs. 3 and 4. Also, the extended coupling modifies 
the stability properties of some trivial phase solutions. For example, the branch (a) is always stable for the 
nearest- neighbor interactions (see top panel of Fig. 1) while it displays an instability window for intermediate 
coupling strengths in the NNN coupling case (see bottom panel of Fig. 1), in a way somewhat reminiscent 
of the findings of [3,4]. 

3. The Variational Approximation 

The Lagrangian of Eq. (4) is 



According to the variational principle, critical points of the Lagrangian (10) correspond to solutions of 
Eq. (4) . This underlies the heuristic justification of the variational approximation ( VA) , whereby an ansatz 
(trial configuration of the wave field) with a finite number of parameters is substituted into the Lagrangian, 
and critical points are then sought for the resulting finite parameter subspace. This approach has long been 
used in various applications, see the review [11] and more recent works, such as Refs. [12-14]. The VA has 
also been used to study nonlocal interactions, see, e.g., Ref. [4], but in the latter case only solutions initiated 
by a single excited site in the AC limit, i.e., solutions of type (a) in terms of Table 1, were studied. To 
the best of our knowledge, the present work for the first time extends the VA to describe not only discrete 
multi-humped solutions with arbitrary phases (with at least three excited sites), but also ones such with 
nontrivial phase distributions. 

3.1. Trivial Phase Distributions 

We start by considering solutions with nontrivial phase distributions. Our objective is to construct ap- 
proximate discrete solitons by means of a real- valued ansatz: 



where S denotes the set of nodes between the first and last excited lattice sites. We define the width W of 
the solution as the number of elements of the set S. Parameters A and _B„ represent the amplitudes, and rj 

^ Generally, the solutions of the nearest-neighbor DNLS may disappear through either pitchfork or saddle-node (in fact, more 
appropriately saddle-center) bifurcations. The configurations considered for this model herein all terminate in the latter type 
of bifurcation, as can be confirmed by comparison with Table 1 of [10] . 




(10) 




(11) 
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is the decay rate. For solutions with odd W (site-centered configurations) the position parameter is hq = 0, 
and for even W (bond-centered configurations) it is no = 0.5. Independent amphtudes _B„ are necessary to 
describe the core of the mult i- humped solutions accurately. 

Since this ansatz is genuinely discrete and includes a purely exponential tail, the corresponding approx- 
imations will only be valid for a small coupling strength, since the solutions become increasingly smooth 
(sech-like) as the continuum limit is approached. On the other hand, solutions with a single excited site in 
the AC limit are only weakly infiuenced by the extended coupling. For this reason, we aim to apply the VA to 
the simplest solution type that "feels" the extended coupling, which corresponds to = 3. The substitution 
of this ansatz into the Lagrangian and calculation of the ensuing sums yields the effective Lagrangian 



(12) 
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The Eulcr-Lagrange equations derived from this Lagrangian are 
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(13) 
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where pi e {A, _B„} are parameters of the ansatz. The equation corresponding to varying amplitude A yields 



2A'^Eo,2r, - 2AEo^r, + e2A ^ -B™^„fc,„ + e-3''e(B_i + Bi)k2 

-f2e-2"e(2Bofc2 -f- (B_i + Bi)ki) = 0, 
and the variation of the Bj's yields, 

B_i{B\ - 1) + e [ Ae-^^'k^ + Ae'^^'ki + ^ 'fc|,„|B|,„|_i J = 0, 
Bo{B^ - 1) + e {2Ae'^'''k2 + k^Bo + fci(B_i + Bi)) = 0, 
B^{Bl - 1) + 6 ( Ae-3"fc2 + Ae-^^k^ + ^ ' k\„,\B^_\^\ J = 0, 

\ mGN / 



(15) 

(16) 
(17) 
(18) 



where the prime over the sum indicates that the m = entry is to be doubled. 

Rather than performing the variation with respect to decay constant 77, we replace it by 77 = In A, where 
A is the corresponding multiplier determined by the dynamical reduction. Furthermore, we set ng = 
rather than using it as a variational parameter. This excludes asymmetric solutions, hence we can also set 
Bi = further reducing the number of parameters. 

Equations (16)-(18) pertain to the localized pattern with W — i, but they can be readily extended to 
other cases. Using these equations, we can approximate all the solutions with trivial phases from Fig. 1, 
i.e., branches (a)-(e), including the single- and double-site solutions, see the top panel of Fig. 2. Indeed, 
the predictions are so good that one cannot see the differences when the bifurcation diagram based on 
the VA is juxtaposed with Fig. 1. A zoom of the saddle-node bifurcation between solutions (c) and (d) is 
shown in the bottom panel of Fig. 2 where the differences arc visible. One might expect the (c) and (d) 
solutions to exchange stability, which would be the case in a standard saddle-node bifurcation in a low- 
dimensional model. However, due to the higher dimensionality of the DNLS model considered here, there 
will be additional eigenvalues that could affect the stability character of the solution. Indeed, for the (c) 
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Fig. 2. Top panels: Numerical solutions of Eq. (1) with trivial phase distributions (solid lines) and the variational solutions 
based on ansatz (11) continued to e = 0.1 (markers). Labels (a)-(d) correspond to those in Table 1. The error \\<j) — at 
these parameter values are (a) 7.4 X 10~^ , (b) .0009, (c) .002, and (d) .001. Middle panel: Zoom of the saddle-node bifurcation 
involving branches (c) and (d) from the bottom panel of Fig. 1. The variational solution (the black dash-dotted line) predicts 
the bifurcation at a slightly smaller value of the coupling strength than the numerical solution (the red dashed line). Both 
solutions are unstable in this parameter region. Bottom panel: Plot of two pairs of isolated eigenvalues versus the arc length of 
the P(e) curve of the middle panel. Here, zero arc length corresponds to the (c) solution at e = and the end of the arc length 
curve corresponds to the (d) solution at e = 0. The eigenvalue pair (I) emerges from the edge of the continuous spectrum, i.e. 

= —1. The eigenvalue pair (II) is responsible for the instability since it has positive real part (i.e. > ) in the entire 
region plotted. 

and (d) solutions there are two pairs of isolated eigenvalues (and one pair at the origin). One of these pairs 
moves along the imaginary axis and becomes real (see curve (I) of the bottom panel of Fig. 2) for some 
critical value of e. However, due to the existence of the other pair of eigenvalues (see curve (11)), which has 
non-zero real part in the parameter region shown, both solutions arc unstable. 

It is obvious that ansatz (11) can only capture solutions with the trivial phase structure, as it is real- 



7 



valued. Nonetheless, it is informative to inspect the AC limit in the framework of the variational equations 
to see what types of solutions are candidates to be approximated. With e = the equations reduce to A = 
and Bj{Bj — 1) = such that the corresponding VA solutions coincide exactly with solutions of the full 
problem (4) with the phases given by either or ±7r, as expected. The fact that the VA and full solutions 
match at e = follows from the choice of the ansatz. This is not the case if the standard ansatz, based on 
the exponential cusp with the single central point is used, cf. Ref. [16]. 



3.2. Nontrivial Phase Distributions 




Fig. 3. (Color online) Left: Phases of the nontrivial solution (f) for e = 0.1 and the same parameter values as in Table 1. Middle: 
Phases at the n = ±1 and n = sites for versus e. Solution (e) corresponds to S_i = d\ = and 6o = i" (red dashed lines) 
while solution (f) corresponds 6o = tt and nontrivial d±i (blue solid curves). The solutions (c) and (f) collide at the critical 
point, e ^ 0.25, where the latter terminates through the ensuing subcritical pitchfork bifurcation. Right: The two eigenvalues 
with the largest real part, found from a solution of Eq. (9) corresponding to solution (e). The eigenvalue corresponding to the 
bifurcation shown in the middle panel is labeled I. However, due to the emergence of a second real eigenvalue pair (II), the 
solution (e) is unstable in all of the parameter region shown. The corresponding spectral picture for the solution (f) is similar, 
but without the curve I. Thus, solution (f) becomes unstable at e fa 0.098. 



In order to capture solutions with nontrivial phases, phase must be added to the ansatz. To motivate our 
choice, we look closer at the solutions with nontrivial phases. Similarly to the trivial-phase ansatz (11), we 
separate the initially excited sites from the others. The outer part of the solution should have exponential 
decay, but with a varying phase. In the left panel of Fig. 3 the phases of the numerically found nontrivial 
solution (f) for e = 0.1 arc plotted against the lattice coordinate. A precise description of each phase would 
make the ansatz too complex, leading to intractable sums in the effective Lagrangian. However, it is clear 
that the solution will be an exponentially localized one, and a coarse approximation for the phases may be 
sufficient. Therefore, motivated by Fig. 3, we assume a constant phase, for n < and another constant 
phase, Ki, for n > 0. For the core of the solution, one might introduce the number of additional phase 
parameters equal to W , one per each site, but, aiming to keep the number of parameters reasonably low, 
we make the following observations based on the dependence of the phases for the core sites as a function 
of e (see middle panel of Fig. 3). We have observed that the phase at the central site (n = 0) stays almost 
constant as the coupling is turned on and the difference between the other points in the core (n = ±1) are 
0±i{e) « O'^i ± 6, where 6*° is the phase at the n-th site in the AC limit, and 6 is a constant which depends 
on parameters of the system (see the middle panel of Fig. 3). Therefore, for the nontrivial-phase solution 
with width = 3 we adopt the following ansatz. 
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ylcxp(iK^i) exp(— ry|n|) for n < — 1, 

Bexp(i((/)° 1 - 6)) for n = -1, 

i'n { Cexp(i^g) for n = 0, (19) 

Bexp(i(0? + 6)) for 71=1, 

Acxp(iKi) cxp(— r/|ri|) for n > 1, 

where real parameters A,B,C represent the amphtude, and the phases are represented by k±i and b. As 
before, 77 is determined via the dynamical reduction. We also present in the right panel of Fig. 3 the instability 
eigenvalues for the trivial phase solution (e). As can be seen from this panel, if it was not for the eigenvalue 
pair (II), this solution would gain stability after collision with the nontrivial phase solution (f) for e > 0.25. 
In particular, this eigenvalue plot in conjunction with the middle panel of the figure illustrating the phases 
of the different branches clearly showcases the existence of a subcritical pitchfork bifurcation, which leads 
to the termination of the nontrivial phase branch (f ) . 

The effective Lagrangian corresponding to the nontrivial-phasc ansatz is the same as in the case of the 
solution with the trivial phases, with the exception that ip appearing in Eq. (14) will be the one of Eq. (19) 
and with B±i = B and Bq = C. 

The variation of the effective Lagrangian with respect to the parameters A, B,C, K-i, ki and b yields, 
respectively. 







)+cos(Ki-0[J)) 



: -2AEo,r, + Beikier^"^ + fc2e~3'')(cos(r") + cos(r+)) + ee"2^Cfc2(cos(K_i - 

+€2A Em^nkrr. + 2A^EQ,2r,. (20) 

rnGS 

: eBk2 cos(6l° 1 + 26 - e^) + e(Cfci(cos(.s-) + cos(.s+)) + A(cos(r-) + cos(r+))(fcie-2'' + fcae-^")) 
+B(S2-l) + eSfco, (21) 
: e(Cfco + Sfci(cos(s^) + cos(s+)) + Afc2e-2''(cos(6l° + k„i) + cos(-6l° + kj))) + C(C^ - 1), (22) 



: eA(Ck2e-^'" sin(K_i - 6*°) + B sin(r-)(fcie-2'? + fcae^^")), 
: eA(Cfc2e~^'' sin(Ki - 6*°) + B sin(r+)(fcie-2»' + fese^^")), 
: eB(2A(sin(r+) - sin(r"))(fcie"2r, _^ fc^g^^'') - Ifca-B sin{el 
— ei?Cfci (sin(s^) — sin(s~)). 



2b -e\)) 



(23) 
(24) 

(25) 



where = k_i - 6 - e\, r+ = ni+b- 6^, = -6^ - 6 + and s+ = -6*? + 6 + 6*^. These equations 
reduce to those for the trivial-phase solutions, displayed in the previous section for k±i = and b = 0^. To 
look at the bifurcation between the trivial- and nontrivial-phase solutions (e) and (f) in another way, we 
can vary the outer coupling parameter k2 in Eq. (5), while keeping e fixed. As mentioned before, it is easier 
to identify the bifurcation through the comparison of phases, therefore we consider the phase difference 
= 6q — Oi^ which has the advantage of being independent of any constant phase (since the solutions 
are gauge invariant). The agreement between the numerically exact solution and the one based on the VA 
is quite remarkable, see the bottom panel of Fig. 4. In this panel, the mirror symmetric nontrivial phase 
solution (i.e., the one with opposite relative phases) has been omitted. 

Using these variational equations, it is possible to approximate the bifurcation point in the bottom panel 
of Fig. 4 which connects solutions (e) and (f) without actually solving Eqs. (20)-(25). Making use of the 
approximations ki = — k_i and k_i = {b -\- 0_i)/2 and Taylor-expanding Eq. (23) leads to 



I B(fcie-2'? + ^26-3';) - 2Cfc2e-2') 
B{kie-^^ + A:2e-3';)/12 - 2Ck2e'^'^i /96' 



(26) 
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Fig. 4. Top panels: The real part (top row) and imaginary part (middle row) of the trivial-phase solution (e) and the one with 
the nontrivial phase (f). The numerically exact solutions (solid lines) are close to their variational counterparts (markers) based 
on ansatz (19), with errors \\<p — il>\\i2 being .004 (e) and .004 (f). The parameters are again the same as in Table 1 and e = 1. 
Bottom panel: Phase difference A8 = 6o — 6i for the numerical (red dashed curve) and variational (black dash-dotted curve) 
solutions for branches (c) and (f) versus the outer-coupling parameter k2 for e = 0.1. The bifurcation here is of the pitchfork 
type, with the mirror-symmetric partner of (f) connecting with solution (c) as well. 



or K-i = 0. Using values of B and C obtained from the (much simpler) equations (15)-(18), we thus produce 
an accurate prediction of the bifurcation, avoiding the need to solve the full system of Eqs. (20)-(25), since 
all terms in Eq. (26) are known. The right hand side of equation (26) vanishes at ^2 ~ 0.538, which is close 
the actual bifurcation value of fc2 « 0.551 in Fig. 4. This is an improvement in comparison to the prediction 
based on Eq. (6), which is k2 ~ 0.5. The deviation from the latter simple leading order prediction can be 
justified by the use of e = 1, whereas the derivation of [2] was valid for small e. 

As a final example, we will briefly consider a four-site solution {W — 4). The corresponding even coun- 
terpart of ansatz (19) is 
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2 



(27) 



^ exp(iKi) exp(— T^l?! + no|) for n > 2, 

where Uq = 0.5. The effective Lagrangian corresponding to ansatz (27) is the same as the trivial-phase 
one, with the exception that once again the ip appearing in Eq. (14) will be of the form of Eq. (27), with 

B_i = B2 = B, Ba = Bi = C, and 

e-'7(2+b-|+2no) _|_ g-»7(4+|j|+2no) 



a2ri _ 1 



(28) 



The extra term e~'''^^+l-'l+^"''^ appearing in the expression for Ej_,^ for the solutions with the even width is 
due to the fact the summation in the Lagrangian is no longer symmetric about the zero site. The variations 
with respect to the parameters A, B, C, ki, & and c yield, respectively, 



= A'Eo.2n - AEo,^ + I] E^^^Kn + eB (fcie-^f^+no) + k2e-'^^'+''°'^j cos{r^) (29) 

rnGS 

+eB (fcie-"(3+"«) + ^26-"^'^+"°)) cos(r+) + eCfc2e-"(^+"°^ cos(r-) + eCfc2e-"'3+"''' cos(r+), 

= 2B{B^ - 1) + 2eBko + eCfci(cos(s^) + cos(s+)) + eCfc2(cos(s-) + cos(s+)) (30) 
+eA (cos(r-)fcie-''(2+"n) +cos(r-)fc2e~'''3+"°' +cos(r+)fcie-''(3+"o) +cos(r+)A^^ 



= 2C(C^ - 1) + eB{ki cos{s~) + ki cos(s+) + ^2 cos(s^ ) + ^2 cos(s+)) + 2eCfco 

+2eCfci cos(6lo + 2c - 6*1) + eAfc2(e-''(2+«o) cos(6'o - k-i + c) + e-''(3+"o) cos(6li - ki - c)), 

= eA{Ck2e~'^^^+"°^^ sin(K„i - 9° - c) + B sm{r-){kie-'^^^+"''^ + fc2e"''(^+"°^), 

= eA{Ck2e''^^^+"°'> sin(Ki - 0i + c) + B sm{r+){kie-'^^^+"°'> + k2e''^^^+"°'>), 

= — ei?C(fc2 sin(s^) + k2 sin(s^) + fci sin(s+) + fci sin(s~)) 

+eBA (sin(r-)fc2e""'^+"") - sin(r+)fc2e""''^+"«' - sin(r+)fcie-"(3+"«) + sin(r-)fcie-"(2+"«)^ 

= ~eCB{k2 sin(s^) + fc2 sin(s+) - fci sin(s+) - fci sin(s^)) - 2fcieC^ sin(6'o - 2c - 6*1) 



(31) 

(32) 
(33) 

(34) 

(35) 



+eAC (sin(r-)fc2e"''(^+"°' -sin(r+)fc2e-''(3+"«)) , 



where ' 



= K_i — 6 — 0_i, TjI" = Ki + & — ^?2, ''^ = «-i — c — 9q and r+ = ki + c — 6*1, = 6'_i —0i+c + b, 
92+ c + b, s~ = 9i — 02 — c + b, s+ = 9^1 — 9q — c + b. Although it is not the goal of this work to 
provide for a list of every possible bifurcation scenario, it is worth mentioning that nontrivial-phase solutions 
connect various types of phase-trivial ones. For example, as depicted in Fig. 5, the nontrivial solution (h) 
connects the phase-trivial ones (g) and (i), with fc2 treated as the bifurcation parameter. Solution (g), with 
phases differences {A9i, A92} = {— tt, — tt}, bifurcates at fc2 ~ 0.345, where the nontrivial-phase solution (h) 
emerges. The phase differences of solution (h) change as ^2 varies. At ^2 ~ 0.894, solution (h) collides with 
and is annihilated by the trivial-phase solution (i), which has phases differences {A^i, A02} = {0, — tt}. I.e., 
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the phenomenology involves a supercritical pitchfork (in the bifurcation parameter fc2) for the emergence of 
the nontrivial phase branch (h) from (g) and a subcritical pitchfork which results from the collision of (g) 
with (i) . This is similar to how asymmetric solutions connect solutions of varying width in DNLS equations 
with higher-order nonlinearities [14,17-19] with the exception that the overall (in)stability of the trivial- 
phase solutions is not affected by collision with nontrivial-phase solutions ^ . The nontrivial-phase solution 
itself undergoes stability change, as can be seen in the bottom panel of Fig. 5. The VA captures this scenario 
remarkably well. Indeed, in the top subpanel of the bottom panels in Fig. 5, the difference between the 
numerical and variational solutions cannot be spotted. Actually, a very strong zoom around the bifurcation 
point is needed to depict the difference (see bottom subpanel) and the relevant error in the identification of 
the critical point is less than 0.4 %. 



4. Conclusions 

We have revisited discrete nonlinear Schrodinger (DNLS) models with extended linear couplings in the 
lattice, developing a new version of the variational approximation (VA) for the models of this type. Using an 
ansatz that coincides with exact solutions in the anti-continuum limit, we were able to accurately describe, for 
the first time, multi-humped solutions, including those with nontrivial phase structures, and the bifurcations 
linking different species of the discrete solitons. Pitchfork bifurcations connecting the solutions with the 
nontrivial and trivial phase structures were identified, which the VA was successful in predicting, thus 
demonstrating its reliability. In particular, the strength of VA analysis is that it provides simple formulas 
to predict such bifurcation points. Like in other settings where the VA is used, a precise evaluation of 
the validity of this approximation is an open question, therefore we relied on the direct comparison with 
numerical solutions of the full problem. It was found that the accuracy of the VA is quite good for small 
lattice coupling strength. 

There remain several other open problems concerning next-nearest-neighbor DNLS equations, such as how 
to compute the stable manifolds in terms of the dynamical reduction, or what the appropriate continuum- 
limit counterparts of these models are and what modifications to the existence and spectral properties of the 
solitary waves the long-range kernels may induce more generally. As concerns the VA, its time-dependent 
version can be also used to predict the linear stability of the solutions, see Refs. [11] and [14], which may be 
another direction for the development of the analysis reported in this work. Finally, bearing in mind some 
of the important consequences of long-range interactions in higher dimensions, such as the stabilization of 
unstable vortices among others [20], consideration of such questions in the discrete realm is a theme of 
particular interest in its own right. 
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^ Nevertheless, along the particular eigendirection of the bifurcation, stability is exchanged between these solutions, as is 
mandated by the pitchfork character of the bifurcation. 
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Fig. 5. (Color online) Top six panels: Profiles of the four-site trivial-phase solutions (g) and (i), which are connected through 
the nontrivial phase solution (h). Bottom two panels: Plots of the phases differences A9i = d\ — do and A62 = 62 — d\ of 
solutions (g)-(i)- At the resolution of the top subpanel the difference between the variational and numerical curves cannot be 
seen. Solution (h) undergoes two stability changes in the parameter region shown here (dashed red and blue solid lines show 
the unstable and stable families, respectively). Bottom subpanel: Zoom of the circled area showing the discrepancy between 
the variational approximation (black dash-dotted lines) and numerically generated curves (red dashed lines). 
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